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Abstract 

We present a density- functional theory based molecular-dynamics study of the struc- 
tural, dynamical, and electronic properties of liquid methanol under ambient condi- 
tions. The calculated radial distribution functions involving the oxygen and hydroxyl 
hydrogen show a pronounced hydrogen bonding and compare well with recent neu- 
tron diffraction data, except for an underestimate of the oxygen-oxygen correlation. 
We observe that, in line with infrared spectroscopic data, the hydroxyl stretching 
mode is significantly red-shifted in the liquid. A substantial enhancement of the 
dipole moment is accompanied by significant fluctuations due to thermal motion. 
Our results provide valuable data for improvement of empirical potentials. 



1. Introduction 



Liquid methanol is of fundamental interest in natural sciences and of signifi- 
cant importance in technical and industrial applications. The liquid phase of 
the simplest alcohol is widely studied, both experimentally and theoretically. 
Among the alcohols, methanol is the closest analog to water. The character- 
istic hydroxyl group allows methanol to form hydrogen bonds that dominate 
the structural and dynamical behavior of the liquid phase. The methyl group 
does not participate in the hydrogen bonding and constitutes the distinction 
with water. This difference is apparent in the microscopic structure of the liq- 
uid, with water having a tetrahedral-like coordination, whereas for methanol 
experiments and molecular simulation suggest a local structure consisting of 
chains, rings, or small clusters. 

The precise quantification of the microscopic structural and dynamical pic- 
ture of liquid methanol has been a long-time subject in both experimental 
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and molecular simulation studies. Recently, a series of state-of-the-art stud- 
ies have been reported. Among these are the neutron diffraction (ND) ex- 
periments of Refs. [1,2]. Simulation studies include work based on empirical 
force-fields [3, 4], mixed empirical and ab-initio interactions [5, 6], and a full ab 
initio molecular dynamics study [7]. The recent ND experiments have provided 
a detailed microscopic picture of the structure of liquid methanol, including 
the pair distribution functions among all atoms. Yet, some of these atom- 
atom distribution functions are still subject to some uncertainty as they are 
obtained indirectly. 

Molecular simulation provides a complementary approach to study the mi- 
croscopic behavior of liquids. Most molecular simulations studies of liquid 
methanol are based on empirical force fields potentials that are designed to 
reproduce a selection of experimental data. Obviously, molecular simulations 
based on these potentials do not provide a picture completely independent 
from experiment. Moreover, the reliability of the results at conditions that are 
significantly different from those where the potential was designed for, may 
be questionable. Density functional theory (DFT) based molecular dynamics 
(MD) simulation, such as the Car-Parrinello molecular dynamics method[8], 
where the interactions are calculated by accurate electronic structure calcu- 
lations provides a route to overcome these limitations. This has been demon- 
strated in studies of liquid water[9,10,ll] and aqueous solvation[12,13,14]. Im- 
portant advantages of DFT-MD over force-field MD are that it intrinsically 
incorporates polarization, that it accounts for the intra-molecular motion and 
therefore allows for a direct comparison with spectroscopy of intra-molecular 
vibrations, and that it yields detailed information on the electronic properties, 
such as the energy levels of electronic states and the charge distribution. In a 
broader chemical perspective it is important to note that DFT-MD is capable 
of the study of chemical reactions in solution, where force-field MD would fail 
completely as it cannot account for the change in chemical bonding. 

Here, we report a DFT-MD study of liquid methanol that addresses the liquid 
structure, the inter- and intra-molecular dynamics, and the electronic charge 
distribution. 



2. Methods and Validation 

Electronic structure calculations are performed using the Kohn-Sham formu- 
lation of DFT. We employed the gradient-corrected BLYP functional[15,16]. 
The choice for the BLYP functional was guided by its good description of the 
structure and dynamics of water [10] where hydrogen bonds are, as in liquid 
methanol, the dominant interactions. Furthermore, it has been shown that 
DFT-BLYP gives a proper description of solvation of methanol in water [13]. 
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The DFT-based MD simulations are performed with the Car-Parrinello method 
[8,17] using the CPMD package[18]. Semi-local norm-conserving Mart ins- Troullier 
pseudopotentials[19] are used to restrict the number of electronic states to 
those of the valence electrons. The pseudopotential cut-off radii were taken 
0.50, 1.11 and 1.23 a.u, for H, O, and C, respectively. The electronic states 
are expanded in a plane-wave basis with a cut-off of 70 Ry yielding energies 
and geometries converged within 0.01 A and 1 kJ/mol, respectively. Vibra- 
tional frequencies are converged within 1 %, except for C-0 and O-H stretch 
modes that are underestimated by 3 % and 5 % compared to the basis-set 
limit values [13]. 

To validate our computational approach we compared results for the gas-phase 
monomer and hydrogen-bonded dimer against state-of-the-art atomic-orbital 
DFT calculations obtained with ADF[20],Q and against the B3LYP and MP2 
calculations of Ref. [21]. The CPMD-BLYP calculations were performed using 
a cubic box with an edge of 12.9 A, with the interactions among the periodic 
images eliminated by a screening technique[17]. Results for geometry and com- 
plexation energy of the dimer are given in Fig. 1 and Table 1. Deviations among 
CPMD and ADF are 1 kJ/mol for the complexation energy, smaller than 
0.01 A for the intra-molecular bonds, and smaller than 0.02 A for the hydro- 
gen bond. This indicates a state-of-the-art accuracy for the numerical methods 
employed in CPMD. Compared to the MP2 and B3LYP results, the BLYP 
bond lengths are slightly longer, with deviations up to 0.03 and 0.05 A for 
the intra- and inter-molecular bond lengths, respectively. Differences among 
BLYP, B3LYP and MP2 complexation energies are within acceptable limits, 
with the BLYP energies smaller by 2-4 kJ/mole. The deviations are similar 
to the comparison between BLYP [10] and MP2Qfor the water dimer and the 
water-methanol dimer [23, 13]. We obtained a zero-Kelvin association enthalpy 
AH°(0) of 10.6 kJ/mol using the B3LYP zero-point energy of Ref. [21]. This 
is in reasonable agreement with the experimental value of 13.2(4) kJ/mol. 
The calculated hydrogen bond length (ro...o = 2.94 A, rn...o = 1-95 A) is 
in good agreement with the experimental values of ro...o — 2.98(2)A and 
r H ...o = 1-96(2) A of Refs. [24] and [25], respectively. 

Current gradient corrected functionals such as BLYP do not account for dis- 
persion forces. For methanol this could be important as attraction to the 
methyl group is fully due to the dispersion force. To estimate the effect of the 
absence of the dispersion we computed the BLYP binding energy of two dimer 
configurations that are sensitive to this: one with methyl groups approaching 
(M-M) and the other with the methyl and hydroxyl group approaching each 



Kohn-Sham orbitals are expanded in an even-tempered, all-electron Slater type 
basis set augmented with 2p and 3d polarization functions for H and 3d and 4f 
polarization functions for C and O. 
2 MP2 limit estimate. See for example [22] 
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other (M-OH). State-of-the art MP2 calculations [26], that incorporate to a 
good approximation the dispersion force, serves as a reference. The dimer 
configurations were taken from Ref. [26] and chosen such that the carbon- 
carbon (M-M dimer) and carbon-oxygen (M-OH dimer) distances were close 
to the peak position of their atom-atom distribution function in the liquid 
state.0 The comparison yields for the M-M dimer values of 2.3 kJ/mol and 
—2.0 kJ/mol for BLYP and MP2, respectively. For the M-OH dimer these 
values are —1.2 kJ/mol and —4.9 kJ/mol, respectively. The too repulsive na- 
ture of the BLYP interaction is consistent with DFT calculations of dispersion 
dominated systems[27,28]. However, although by far not insignificant, the mag- 
nitude of the deviation is much smaller than the hydrogen-bond interaction 
and of the same order of magnitude as the error in the latter. It can therefore 
be argued that for a study of liquid methanol on the accuracy level of BLYP, 
neglecting the dispersion interaction is acceptable. 

In Ref. [13] we have shown that for the gas-phase monomer CPMD-BLYP 
vibrational frequencies are in excellent agreement with ADF results and, com- 
pared to experiment, underestimate allmost all modes by ~ 10%, a known 
feature of the BLYP functional. 

Overall, we conclude that our level theory is satisfactory in comparison with 
experimental and other theoretical gas-phase data. 



3. Liquid 

Liquid methanol was modeled by 32 molecules in a periodic cubic box with 
an edge of 12.9 A, reproducing the experimental density of 0.791 g/cm 3 at 
293 K[29]. The temperature was fixed at 293 K using the Nose- Hoover ther- 
mostat [30]. The fictitious mass associated with the plane- wave coefficients is 
chosen at 900 a.u., which allowed for a time step in the numerical integration of 
the equations-of-motion of 0.145 fs. The system was equilibrated for 1 ps from 
an initial configuration obtained from a force-field simulation. Subsequently, 
we gathered statistical averages from a 6 ps trajectory. 

3.1 Structure 

In Fig. 2 we have plotted the most characteristic atom-atom radial distribution 
functions (RDFs), i.e. the hydrogen bonding O-O, O-Hq, and Hq-Hq RDF and 

3 Configuration from Ref. [26]. M-M dimer: geometry M with rcc = 3.75 A. M-OH 
dimer: geometry I with r^o = 3.50 A. Both monomers in the dimer are kept fixed 
to their isolated geometries. 
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the C-0 RDF. For comparison we also plotted results of recent ND results [1] 
and the peak positions obtained using Haugney's empirical potential [3] , the 
latter being considered one of the most accurate empirical force fields to date. 

The pronounced structure in the first three RDFs are a clear indication of the 
presence of hydrogen bonds. Comparison with the experimental data shows 
that the positions of the first peaks match within the statistical error for the 
0-0 and 0-H o RDFs and is slightly smaller (« 0.1 A) for the H -H RDF. 
The height of the first peak is in good agreement for the O-Hq and Hq-Hq 
RDFs, that both can be determined accurately from the ND data. However, 
the 0-0 RDF shows a calculated first peak height that is significantly lower 
than the experimental result. Given the small system of 32 molecules in our 
simulation, the discrepancy could well be a system-size effect. On the other 
hand, the indirect way by which the the 0-0 RDF is extracted from ND data 
could yield an overestimate of the 0-0 correlation. Comparison with force 
field results [3], that yield significant higher peak values for the 0-0 and O- 
Ho RDFs, suggests that the Haugney potential overestimates the hydrogen 
bonding structure in the liquid, in line with the observation of Ref. [1]. 

The number of H-bonds as calculated by integrating the 0-Ho and 0-0 RDFs 
up to the first minimum, and using the geometrical criterion of Ref. [3] , yields 
values of 1.9, 2.0, and 1.6, respectively. This is in good agreement with the 
experimental ND results of Ref. [1] yielding 1.8 and 1.9 obtained by integrating 
the 0-Ho and 0-0 RDFs. Applying the geometrical criterion to the Haugney 
force-field simulation[3] yields a slightly higher value of 1.9. 

The hydrogen bonding in the liquid phase is accompanied by an elongation of 
the OHq bond of 0.15 A. A direct comparison with the experimental results 
for this change in the geometry of the methanol molecule is rather difficult 
due to the large spread in the reported values. However, similar change in the 
geometry is observed in the DFT-MD study of liquid methanol reported in 
Ref. [7] and in ab initio studies of small methanol clusters[21]. 

The calculated C-0 RDF is in reasonable agreement with the ND results, with 
the overall shape well reproduced but the first peak clearly less pronounced 
than in the ND result. This is consistent with the, in the previous section 
found, underestimation of the BLYP binding energy of the M-0 dimer with 
a C-0 distance at the RDF peak position. This is due to the absence of 
the dispersion interaction in BLYP. However, the absence of the dispersion 
interaction clearly does not lead to a completely distorted C-0 positional 
correlation. The comparison of the calculated C-C RDF with the ND result 
(not plotted) is very similar. 
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3.2 Dynamics 



The time scale of the present simulation (6 ps) allows for an analysis of the 
short-time dynamics of liquid methanol. Figure 3 shows the power spectrum of 
the velocity auto correlation function (VACF) of the hydroxyl hydrogen. For 
comparison we have also plotted the calculated 200 K monomer spectrum of 
Ref. [13]. The three distinct peaks correspond to the OH stretch (3100 cm" 1 ), 
C-O-H bend (1600 cm" 1 ), and the CO stretch (1000 cm" 1 ). The broad feature 
below 1000 cm -1 indicates the librational-translational (500 cm" 1 ) modes of 
the methanol molecules. Compared to the gas phase, the liquid OH stretch 
mode has red-shifted by approximately 200 cm -1 and broadened consider- 
ably. On the other hand, the C-O-H bending is blue-shifted by approximately 
70 cm" 1 . The observed shifts and broadening are characteristic for hydro- 
gen bonded liquids and also observed in the spectrum of water or hydrated 
methanol. The calculated shifts compare reasonably well with experimental 
infrared spectra[31] that yield values of —354 cm" 1 and +78 cm" 1 for the O- 
H stretch and C-O-H bend. The calculated positions and shifts of the modes 
match within statistical errors with those of methanol in aqueous solution[13] 
determined using the same computational approach. This indicates that the 
intra-molecular dynamics of methanol is affected in a similar way by an aque- 
ous environment and a methanol environment. 

The diffusion constant D is a key measure of the collective dynamics. In view 
of the limited length of the calculated trajectory we can only provide a rough 
estimate. From the mean square displacement of the oxygen atoms we obtained 
D = 2.0 ± 0.6 x 10" 9 m 2 /s, in reasonable agreement with the experimental 
value of 2.42 ± 0.05 x 10~ 9 m 2 /s [32]. 



3. 3 Electronic properties 

As the electronic structure is an intrinsic part of a CPMD simulation, detailed 
information on the electronic charge distribution is obtained. To quantify the 
charge distribution we used the method of maximally localized Wannier func- 
tions that transforms the Kohn-Sham orbitals into Wannier functions whose 
centers (WFC) can be assigned with a chemical meaning such as being asso- 
ciated with an electron bonding- or lone-pair (LP) [33]. 

We calculated the positions of the WFCs for the monomer, the dimer, and 6 
independent configurations of the liquid simulation. Table 2 lists the (average) 
distances of the WFCs associated with the oxygen electrons. Most notably is 
the small but significant shift of 0.024 A for the OH bond WFC towards the 
oxygen atom when going from the monomer to the liquid. At the same time, 
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one of the LP WFCs shifts away from the oxygen by 0.023 A. These changes 
should be considered a manifestation of the hydrogen bonding and the induced 
polarization among the dipolar methanol molecules in the liquid state. 

To quantify the change in the charge distribution in a single number we cal- 
culated the molecular dipole moment assuming the electronic charge to be 
distributed as point charges located on the WFCs. For liquid water it has 
been shown that such a partitioning of the charges over the molecules yields 
a unique assignment of the WFCs over distinct moleculesfll]. From Table 3, 
that lists the values for the monomer, dimer, and liquid, we observe a sig- 
nificant enhancement of the dipole moment going from the monomer via the 
dimer to the liquid. A comparable liquid-state value of 2.39 D has been ob- 
served in a coupled empirical and ab initio MP2 study [5]. Note that the value 
of the dipole moment is somewhat larger than in the Haugney[3] (2.33 D) or 
AMBER[34] (2.2 D) force field. A second important feature of the electronic 
charge distribution in the liquid is its fluctuating character due to the ther- 
mally driven configurational changes. In Fig. 4 we have plotted the calculated 
distribution of the dipole moments in the liquid phase. It shows that there is 
a significant variation ranging from 1.7 D to 3.5 D. 



Conclusions 

We have demonstrated that ab initio MD is a valuable approach to study 
the structural, dynamical, and electronic properties of liquid methanol. The 
calculated pair distribution functions involving the hydroxyl hydrogens corre- 
late well with recent state-of-the-art neutron diffraction experiments of Soper 
and co-workers. It confirms their finding that one of the benchmark empiri- 
cal potentials overestimates the hydrogen bonding structure. The calculated 
oxygen-oxygen radial distribution function shows significantly less structure 
than the experimental neutron diffraction results. Currently we are studying 
a larger simulation sample to see whether this is due to the small system size 
in the present calculation. It could also be a result from an inaccuracy in the 
experimental result that is obtained in an indirect way. Results for the dimer 
binding energies and the oxygen-carbon RDF suggest that the absence of the 
dispersion interaction is notable but has no major impact. Comparing the vi- 
brational spectra of the liquid phase against that of the gas phase monomer 
shows a significant red shift of the O-H stretch accompanied by a smaller blue 
shift of the C-O-H bend mode, in reasonable agreement with experimental 
observations with the O-H shift somewhat underestimated in our calculation. 
We quantified the electronic charge distribution using a Wannier function de- 
composition. A small but measurable shift of the positions of the Wannier 
function centers when going from the gas-phase to the liquid is accompanied 
by a substantial enhancement of the dipole moment. Moreover we have found 
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that in the liquid the dipole moment fluctuates significantly with variations 
up to half the average magnitude. The latter suggest that the assumption 
made in empirical potentials using a fixed dipole moment is a strong simplifi- 
cation. The present results may be considered valuable data for improvement 
of empirical potentials for the study of liquid methanol. 
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Table 1 

Complexation energies (kJ/mol) of methanol dimer shown in Fig. 1. Numbers are 
bare values without zero-point energy corrections and without entropy contribu- 
tions. 



CPMD-BLYP 


ADF-BLYP a 


B3LYP 6 MP2 C 


16.4 


17.3 


20.6 18.4 



a Refs. [20] 

6 B3LYP/6-311+G(3df,2p) method. B3LYP/6-311+G(d,p) optimized geometries. 
From Ref. [23]. 

c G2(MP2) method. MP2(full)/6-311+G(d,p) optimized geometries. /.From 
Ref. [23]. 

Table 2 

Electronic charge distribution in terms of Wannier function centers. rf(LP) denotes 
the average distances between a lone pair WFC and the O nucleus. d(OH) and 
d(OC) denote the (average) distances between the covalent WFC along the O— H 
bond and the O— C bond with the O nucleus, respectively. All distances are given 
in A. Statistical errors for the liquid data are around 0.002. 





d(LP) 


d(LP) 


d(OH) 


d(OC) 


Monomer 


0.305 


0.305 


0.533 


0.562 


Dimer 


0.316 


0.306 


0.522 


0.561 


Liquid 


0.328 


0.309 


0.509 


0.561 



Table 3 

Dipole moment. Experimental value is given in parentheses. Data for the liquid 
phase were obtained by averaging over 6 configurations of the MD simulation. Sta- 
tistical errors are in the order of some units in the last digit. 



M (D) 


Monomer 


1.73 (1.69°) 


Dimer 


2.03 


Liquid 


2.54 



Microwave study, ref. [35]. 
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0.962 



Fig. 1. Optimized geometry of the methanol dimer. Selected distances (A) and 
angles (degrees) are given for three computational methods: CPMD-BLYP (top, 
this work), ADF-BLYP[20] (second, this work) and B3LYP[21] (third). The MP 2 
results of Ref. [21] are within 0.01 A of the B3LYP result. 
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Fig. 2. Calculated hydrogen-bonding and C-0 radial distribution functions (solid 
lines). Dashed line indicate neutron diffraction results of Ref. [1]. Crosses indicate 
position of the first peak of the RDFs obtained by Haugney et al. using an empirical 
force field [3]. 
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Fig. 3. Calculated power spectrum of the VACF of the hydroxyl hydrogen for an 
isolated methanol at T=200 K (solid line, from Ref. [13]) and liquid methanol at 
T=293 K (dashed line). 
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Fig. 4. Distribution of the molecular dipole moment in liquid methanol, obtained 
from 6 independent liquid configurations. 
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